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Abstract. A scheme stemming from the use of pseudospectral approximations to spatial 
derivatives followed by a time integrator based on trigonometric polynomials is proposed 
for the numerical solutions of the coupled nonlinear Klein-Gordon equations. Numerical 
tests on one- and three-coupled Klein-Gordon equations are presented, which are geared 
towards understanding the accuracy and stability, and demonstrating the efficiency and 
high resolution capacity in application. 

1. Introduction 

The characteristics of nonlinear phenomena in various physics fields as, e.g., fluid dynam- 
ics, laser and fiber optics, solid state physics, plasma physics, chemical physics, reaction 
kinematics and etc., can be mathematically described by nonlinear evolution equations pQ. 
Among them of great physical significance are the equations that possess soliton solutions. 
In recent decades a variety of methods, such as the inverse scattering method [1,2J, bilinear 
transformation |12j and etc., have been developed for obtaining the exact solutions of soli- 
ton type equations. In parallel with the analytical treatment a surge of studies have been 
devoted to the numerics of these equations, which is a topic of great importance in applied 
science. In the present work the numerics of coupled nonlinear Klein-Gordon equations 
governing waves in a dispersive media is going to be examined. 

TV-coupled nonlinear Klein-Gordon equations under consider take the general form [3lfTT] : 

(i.i) 
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N 
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with k = 1, 2, . . . , iV, for ^(x, t) and Q(x, t) sufficiently differentiate functions. Note that 
there exists an invariant of (|1.1I) - (|1.2I) : 

(i.3) w = jf |E[w*) 2 + ( a ^*) 2+ ^]-(E^) +\Q 2 } dx ' *^°> 

and it suffices to refer the above conserved quantity as energy. Such conservation law can 
be justified from multiplying both sides of fjl.ll) by 2dt^k an d, respectively, both sides of 
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(II. 2p by Q, integrating over R and then summing up for k = 1, 2, . . . , N. Bilinear form and 
one-soliton solutions of (jl.ll) - (jl.2l) were investigated in [3|ITT]. and the complete integrability 
was also constructed from P-analysis [H[T5]. Along the numerical aspects of (jl.ll) - (jl.21) . to 
our best knowledge there are few results derived in literature. 

The goal of this paper is to propose an efficient and accurate numerical scheme for solving 
(jl.ll) - (jl.21) . The key point in designing the scheme is based on applying Fourier pseudospec- 
tral approximations to the spatial derivatives, followed by applying a time integrator based 
on trigonometric polynomials, i.e., the so-called Gautschi-type or Deuflhard-type integrator 
(see e.g. [8, 10J) in phase space to the temporal discretization. The resulting scheme is fully 
explicit, symmetric in time, spectral-order of accuracy in space and second-order of accu- 
racy in time, easy to implement and with less memory demand. Moreover, numerical tests 
demonstrate that the scheme is stable and conserves the energy defined by (11.31) very well, 
which are two favorable properties desired for a scheme in long-time simulation application. 
Note that the similar technique has been used for standard non-coupled Klein-Gordon equa- 
tions [5j, and for other classes of schemes, as for instance the finite-difference, decomposition 
and etc., for non-coupled Klein-Gordon equations we refer the readers to [l6l[7l[9l [T3l[T4] and 
references given therein. 

The rest is as follows. In Section [2] efficient numerical scheme is proposed. In Section 
[3] numerical results are reported for accuracy and stability tests of the scheme, and its 
application in numerically studying the dynamics of soliton-soliton collisions in one- and 
three-coupled Klein-Gordon equations. Finally some conclusions are drawn in Section [H 

2. Numerical scheme 

In this section we shall propose the efficient numerical discretization. The initial condi- 
tions are assumed to be of the form: 

(2.1) il> k (x,t = 0)=rl>£\x), dtM^t = 0)=^\x), Q(x,t = 0) = Q(°\x), xEi 

In practice we truncate the whole space problem into an interval [a, b] with periodic boundary 
conditions 

i/j k (a, t) = ij) k (b, t), d x ip k (a, t) = d x ^ k {b, t), Q(a, t) = Q(6, t), t > 0. 

Choose mesh size h := (b — a)/M with M an even positive integer, time step r > 0, 
and denote the grid points with coordinates (xj,t n ) = (a + jh, nr) for j = 0, 1, . . . , M and 
n — 0,1,.... Let an d be the approximations of ip k (xj,t n ) and Q(xj,t n ), and 

denote by ij) k and Q n the solution vectors with components('0/ c )J and Q 7 - respectively. 

2.1. Discretization for (11.11) . Assume 

M/2-1 

(2.2) ip k (x, t) w ^ ($k)i(t) exp (im{x - a)) , 

J--M/2 

for a < x < b and t > 0, with ^ = 2tt//(6 - a) (I = -M/2, . . . , M/2 - 1) and (^) z (t) the 
discrete Fourier transformation coefficient of the Ith mode, 

^ M-l 

(2.3) (^ fc )z(t) = M*j,t) exp (-ini(xj - a)) , 1 = -M/2, . . . , M/2 - 1. 

3=0 
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Approximating the spatial derivatives in (11.11) by Fourier pseudospectral discretization [16], 
i.e., 

M/2-1 

-d xx ^ k {x,t) « ^2 vK^k)i(t) exp (ifj,i(x - a)) , 

J--M/2 

and noting orthogonality of the Fourier functions, we obtain the following second-order 
ODEs in phase space, for / = -M/2, . . . , M/2 - 1, 

(2.4) + (m? + 1) («(*) - (AM*) = o, 

with fk(x, t) — 2 (j2p = i ifjp(x, t) + Q(x, t)^j iik(%, t) and (fk)i(t) defined in an analogous way 
as (|2.2I) . The analytical solutions of the above second-order ODEs can be written explicitly 
thanks to variation-of-constants formula. For t n (n=0,l,. . . ) a given time, 

(tpk)i(t) = (tpk)i(tn) cos (\i(t- t n )) + X^ 1 ^(^ k ) l (t n )siri(Xi(t- t n )) 

(2.5) +A^ r(A)iWsin(Ai(t-fi))da 3 



with A/ = y /i^ + 1. 

The approximations of ij)k{x,ti) are achieved from (|2.5I) with n = and £ = ti = r to- 
gether with the initial conditions (|2.1j) and approximating the integrals by standard trape- 
zoidal rule. For n = 1,2,..., adding (|2.5I) with £ = t n+ i = t n + r to its evaluation at 
t — t n -i — t n — r, and then approximating the integrals via trapezoidal rule, we get the 
three-term recurrence: 

fe)i(Wi) = - O0fe)z(*n-i) + 2 (^k)i(t n ) cos (A/r) 

+ V 1 / f(X)|(*n + 5) + (A)z(*n " s)l sin (A«(r - s)) ds 
Jo 1 J 

(2.6) « - fi/S^tn-i) + 2(^(t n ) cos (^ T ) + rX^(K)i(t n ) sin(V). 
2.2. Discretization for (11.21) . Again, assuming 

M/2-1 

Q(x,t)w ^ Qi(t) exp(i^(x - a)) , 

l=-M/2 

for a < x < 6, t > 0, and approximating the spatial derivative in (|1.2I) by Fourier pseu- 
dospectral discretization [IB], 

M/2-1 

(2.7) d x Q{x,t) « ^ mQi{t) exp(z/x/(x - a)) , 

l=-M/2 

we have the following first-order ODEs in phase space, for / = —M/2, . . . , M/2 — 1, 
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with g(x,t) = — 2 J2p=i V^O 3 ^ t)- For t n (n = 0, 1,. . .) a given time, integrating (|2.8p from 
t — t n to t — t n +i = t n + r and approximating the integrals via trapezoidal rule, we get, 

Qz(*n+i) = exp (ifiir) Qi(t n ) + gi(t n+1 ) - exp (i/^r) <#(£ n ) 



rtn+i 

+ z/i; / exp (i^i(t n+1 - t)) gi(t)dt 

Ju. 



(2.9) « exp (i/x z r) Qz(*n) + + x j 5/(*n+i) + ~ X j ex P (w r ) 

2.3. Detailed numerical scheme. Choosing = an( ^ = the 

detailed numerical scheme is as follows: for n = 1, 2, . . . and j = 0, 1, . . . , M 

M/2-1 



(2-10) = E (^W(^)> 

l=-M/2 V 7 

M/2-1 / . 

(2-n) Ql= E (w-pf^ 1 )' 

l=-M/2 V 7 

where, 

(2.12) (^ = (^°^cos(A,r) + (V^)^" 1 sin ( V) + ^ (/£), sin (A z r) , 

(2.13) tyj^), = -Wf\ + 2(^jf) J cos (V) + r\l l (f») l sin(A,r), n = 1, 2, . . . , 
(5^), = exp (i/nr) (Q^ + + l) 5^), 

(2.14) + P^-l) exp (i^r) (F)„ n = 0, 1, . . . . 



Here, the vector /£ = [(/ fc )£ , , . . . , (/ fc )^] T and s» = <ff , . . . , <^] T are defined by 

(N 2 \ JV 2 

E ((W) + q? (^)" , ^ = -2 E ((^), n ) • 
P =i / P =i 

If the energy defined by (|1.3p is of interests, then the approximation of dtipk( x j,tn)i 

M/2-1 . . 

l=-M/2 V 7 

for j = 0, 1, . . . , M and n = 1, 2, . . . can be obtained from 

(2.16) - (yf\ = 2A r 1 (^) / sin(A,r) , 

which is achieved by subtracting (|2.5p with t = t n +i = t n + r from its evaluation at 
£ = — t n — t, and applying trapezoidal rule to the integrals. 

The scheme (|2.10l) - (|2.16l) is fully explicit, symmetric in time by noting that it is unchanged 
if we interchange r <H> — r and n <H> n + 1, and quite efficient in implementation thanks to 
FFT. It is of spectral-order of accuracy in space, i.e., it converges exponentially fast as 
mesh size refined smaller, which is an expected property for the spectral-type discretization. 
Moreover, as shown by numerical experiments reported in the next section, the scheme is 
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Table 1. Accuracy tests results: (i) upper part for discretization error in 
space, under r = 0.0001; (ii) middle part for discretization error in time, 
under h = 1/8 and (iii) lower part for conserved quantity analysis, under 
h = 1/4 and r = 0.02. 





h = 1/2 


h = 1/4 


h = 1/8 


h = 1/16 


e(60) 


1.1677E-1 


2.8638E-6 


1.7098E-6 


1.0244E-8 




r = 0.04 


r = 0.02 


r = 0.01 


r = 0.005 


e(60) 


1.1654E-1 


2.9405E-2 


7.3659E-3 


1.8423E-3 


t 


t = 50 


t = 100 


t = 150 


t = 200 


E(t) 


0.67890052 


0.67890052 


0.67890050 


0.67890051 



of second-order of accuracy in time, stable and conserves the energy defined by (11.31) very 
well. 



3. Numerical results 

Numerical examples are reported in this section towards understanding the accuracy 
and stability of the numerical scheme fl2.10l) - fl2.16l) . and applying it to study soliton-soliton 
collisions in one- and three-coupled nonlinear Klein-Gordon equations. 

3.1. Accuracy tests and stability study. It is known that the Af-coupled nonlinear 
Klein-Gordon equations admit the following analytical one-soliton solutions [3]: 

j\ -\- c ( x — ct \ 

(3.1) ^( c ,a k )(x,t) = a *y sedl vTT^py ' 

(3.2) QcOM) = ~~ ~~t sech 2 X ~ Ct 



i \VT^ 

with \c\ < 1 an arbitrary constant, and coefficients a& satisfying J2k=i a k ~ ^- ^° ^ es ^ ^ ne 
accuracy, we solve the one-coupled nonlinear Klein-Gordon equations (i.e. N = 1 in (11.11) - 

(OJ> ) on [-24, 104] for < t < 200, with initial data ^\x) = *Jj {0A q(x,t = 0), ^\x) = 

9t^(0A, i)( x : t — 0) an d Q^°\x) = Qo.4(x,t = 0). To quantify the numerical results, the 

error function e{t) is defined as the maximum error, e(t n ) := maxj (^i)^ — ^i(xj,t n) 



+ 



maxj Q^Xjit^ 

To test the discretization error in space, we choose a very fine time step r = 0.0001 such 
that the error from time discretization is negligible compared to the spatial discretization 
error. Similarly, we choose a very fine mesh size h = 1/8 to eliminate the spatial discretiza- 
tion error for testing the discretization error in time. Also, mesh size h = 1/4 and time step 
r = 0.02 are chosen to test the conservation of conserved quantity (jl.3p . The results are 
listed in Table [TJ To study the stability of the scheme, the same initial conditions as chosen 
previously are used, to which we add white Gaussian noise with signal-to-noise ratio 50dB. 
The results are depicted in Figure [TJ 

Based on these results, the following are drawn: 
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Figure 1. Comparisons between numerical one-soliton results (solid line) 
from initial data perturbed by adding white Gaussian noise with signal- 
to-noise ratio 50dB and results (dashed line and red colored online) from 
initial data without perturbation, under h = 1/4 and r = 0.02. Wave fronts 
depicted for t = 0, 50, 100, 150 and 200 (from left to right). 



(1) The numerical scheme (12. 10}) - (12. 16f> is of spectral-order of accuracy in space, i.e., its 
convergence order in space is higher than any polynomial order, and of second-order 
of accuracy in time. 

(2) It conserves the invariant defined by (11 .3|) very well (up to seven significant digits in 
the reported example) over a long-time simulation. 

(3) It is numerically stable, by which we mean a minor perturbation in the initial data 
dose not bring in significant difference to the results over a long-time simulation and 
no numerical blow-up occurs. 

3.2. Applications on soliton-soliton collisions. First we apply the method to one- 
coupled nonlinear Klein-Gordon equations, i.e., N = 1 in (ll.ip - (11.2p . with initial date: 

(3.3) -01 O) 0) = -0(0.6, 1) + x 0>t = 0) + 0(_ O .25, i)(x-x ,t = 0), 

(3.4) ^i\x) = <9t^(o.6, i)0 + x ,t = 0) + <9f0(- O .25, 1)0 -x ,t = Q), 

(3.5) Q(°) 0) = Qo.eO + so, * = 0) + Q-0.25O - x , t = 0), 

with a dislocation parameter xo > to measure the separation of two one-soliton waves 
at initial time level. Some results from numerical experiments are shown in Figure [2] and 
Figure [31 These results indicate that when two solitons are of nearly completed separation 
at initial time level, after the collision the solutions remain to propagate in soliton pattern as 
the time proceeds and there is no visible wave structure emitted. On the other hand, when 
two solitons are not quite separated at initial time level, after the collision the emission of 
new waves is conspicuous. 

Next we present an example of applying the method to three-coupled nonlinear Klein- 
Gordon equations, i.e., N = 3 in (ll.ip - (ll.2p . and the initial conditions are taken as 

(3.6) ^i° } 0) = -0(0.6, l/v^O + x ^ = 0) + ^(-0.25, -1/2)0 -x ,t = Q), 

(3.7) 4° } 0) = 0(o.6, -iA/2)0 + = 0), 

(3.8) ^3° } 0) = ^(-0.25, y/S/2)( x - X 0^ = 0) ? 
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Figure 2. Numerical results for soliton-soliton collisions in one-coupled 
Klein-Gordon equations. Left column for solution evolution and right col- 
umn for its contours (color online). Computations are carried out on [—24, 40] 
under h = 1/4 and r = 0.02. Initial conditions are chosen as (|3.3l) - (j3.5l) with 
xq = 8 such that initially two one-soliton waves are well separated. 



(3.9) ^i\x) = 0*^(0.6, 1/V2)( X + ^0,* = 0) + <9^(_ .25, -i/2) (x -x ,t = 0), 

(3.10) ^2 ] ( x ) = ^(0.6, -i/V2)( x + ^o,£ = 0), 

(3.11) ^\x) = Qfc^(_ .25, >/3/2)( X - ^0,* = 0), 

(3.12) Q(°) (x) = Q 0>6 (x + x 0l t = 0) + Q- .2s(x - x ,t = 0), 

again with a dislocation parameter xo > to measure the separation of two one-soliton 
waves at initial time level. Some results from numerical experiments are depicted in Figure 
S] and Figure [5j in which initially two soliton waves are well separated and only a wave front 
moving towards the right x-axis exists in ^2 while only a wave front moving towards the 
left x-axis exists in ^3. It shows that after the collision a new wave front moving towards 
the left x-axis appears in ^2 and a new wave front moving towards the right x-axis appears 
in ^3. Also, a conspicuous change in the amplitude of waves is observed, see Figure E 
Results for two solitons not separated at initial time level are quite similar as in one-coupled 
Klein-Gordon equations. 

The results presented in this section indicate that the dynamics of waves governed by 
coupled nonlinear Klein-Gordon equations is a rather complicated issue. The results also 
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X X 

Figure 3. Contours (color online) of numerical results for soliton-soliton 
collisions in one-coupled Klein-Gordon equations. Computations are carried 
out on [—32, 32] under h = 1/4 and r = 0.02. Initial conditions are chosen as 
(I3.3p - (l3.5p with xq = 1 such that initially two one-soliton waves are almost 
at the same location. 

demonstrate the efficiency and high resolution of the proposed method for numerically 
studying the coupled nonlinear Klein-Gordon equations. 

4. Conclusions 

A numerical scheme, which is based on the application of pseudospectral approximations 
to spatial derivatives followed by a trigonometric integrator to temporal discretization in 
phase space, was proposed for solving coupled nonlinear Klein-Gordon equations. The 
soliton solutions arising from one-coupled Klein-Gordon equations were examined to test 
the accuracy and stability. Also, application results of studying soliton-soliton collisions in 
one- and three-coupled Klein-Gordon equations were reported as examples to demonstrate 
the efficiency and high resolution of the scheme. 
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Figure 5. Numerical results of ipi, ip2 and ^3 for soliton-soliton collisions in 
three-coupled Klein-Gordon equations: solid line corresponds to time t — 200 
and dash-dot line (red colored online) corresponds to time t — 0. Parameters 
are same as Figure [H 
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